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Columnar crystals contain defects in the form of vacancy/interstitial loops or strings of vacancies 
and interstitials bounded by column "heads" and "tails". These defect strings are oriented by the 
columnar lattice and can change size and shape by movement of the ends and forming kinks along 
the length. Hence an analysis in terms of directed living polymers |lj is appropriate to study their 
size and shape distribution, volume fraction, etc. If the entropy of transverse fluctuations overcomes 
the string line tension in the crystalline phase, a string proliferation transition occurs, leading to 
a supersolid phase Q. We estimate the wandering entropy and examine the behaviour in the 
transition regime. We also calculate numerically the line tension of various species of vacancies and 
interstitials in a triangular lattice for power-law potentials as well as for a modified Bessel function 
interaction between columns as occurs in the case of flux lines in type-II superconductors or long 
polyelectrolytes in an ionic solution. We find that the centered interstitial is the lowest energy defect 
for a very wide range of interactions; the symmetric vacancy is preferred only for extremely short 
interaction ranges. 
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I. INTRODUCTION 



The physics of columnar crystals is relevant to the Abrikosov lattice of flux lines in Type-II superconductors and 
liquid crystalline materials like concentrated phases of long polymers or discotics. The stability of the columnar 
crystal has been investigated, and various mechanisms proposed for its melting. Conventional melting, which arises 
when phonon displacements reach a fixed fraction of the lattice constant, can easily be located via the Lindemann 
criterion Melting destroys the two-dimensional crystalline order perpendicular to the columns leading to a 

nematic liquid of lines or columns, which is entangled at sufficiently high densities. 

Crystal defects play an important role above the melting transition. If edge dislocations in the crystal proliferate, 
they drive the shear modulus to zero, leading to a liquid-like shear viscosity. However, dislocations alone cannot 
destroy the six-fold orientational order of the triangular lattice in a two-dimensional cross-section. Thus, provided 
disclination lines do not also proliferate, the resulting liquid of lines is hexatic, not isotropic The screw component 
of the unbound dislocations leads to entanglement. A finite concentration of unbound disclinations superimposed on 
the hexatic liquid leads to isotropic in-plane order. 

Another kind of transition is brought about by vacancy /interstitial line defects in columnar crystals composed of 
long, continuous lines. As discussed in Ref. [H, under suitable conditions (such as high field and small interlayer 
coupling in layered superconductors), it can become favourable for these line defects to proliferate. If this happens 
at a temperature Td below the melting temperature T™, then the phase that exists between Td and Tm will be 
simultaneously crystalline and highly entangled. In the boson analogy of an aligned system of lines, where the lines 
represent two-dimensional bosons traveling in the "time-like" axial (z) direction such a phase is analogous to 
the supersolid phase of the bosonic system which incorporates vacancies and interstitials in its ground state. This 
entangled solid melts into an entangled liquid or an entangled hexatic at even higher temperatures. 

The proliferation of vacancy or interstitial strings could also affect a crystal-to-hexatic transition mediated by 
dislocations. Dislocations in the columnar crystalline geometry are normally constrained to lie in the vertical plane 
formed by their Burger's vector and the z-axis, because a dislocation in a two-dimensional cross-section can move 
along the columnar axis only through glide parallel to its Burger's vector. Transverse motion (climb) would require it 
to absorb or emit vacancies or interstitials. This becomes possible in the supersolid phase, thus allowing dislocation 
loops to take on arbitrary non-planar configurations which would have to be included in the treatment of Ref. Q to 
study melting out of a supersolid phase |^ . 

Vacancy/interstitial strings in a columnar crystal tend to be lines themselves because of the continuity of the 
columns. If the columns are constrained to be continuous across the entire sample (as is the case for vortex lines in 
Type II superconductors), these defects must either thread the entire sample (Fig. ^ or appear in vacancy/interstitial 
pairs forming loops (Fig. ||) ||^. The situation is different, however, for finite-length polymers, or columns of discotic 
liquid crystal molecules which can break and reform freely. As illustrated in Fig. |3|a, a slice through a low temperature 
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configuration in a polymer columnar crystal (with translational order perpendicular to the column axis but not parallel 
to it) would consist of tightly bound polymer "heads and tails" . At higher temperatures, however, the heads and tails 
will separate, either moving apart to form a vacancy string or sliding past each other to form a line of interstitials 
(Fig. gb) 0. In columnar discotic crystals with similar translational order, "heads" and "tails" are absent at low 
temperatures, but appear spontaneously when vacancy and interstitial strings are excited (Fig. (Head and tail 
defects appear superficially like dislocations in the cross sections shown in Figs. ^ and |[ A three-dimensional analysis 
of lines and columns in neighbouring sheets like that shown in Figs. |^ and ^ is necessary to clearly reveal that these 
are strings of vacancies and interstitials.) 
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FIG. 1. Vacancy string rd{z) (thick dashed curve) mean- i i i i i i i i i 

dering through a columnar crystal. Dashed lines represent FIG. 2. Vacancy-interstitial loop in a columnar crystal, 

columns just above or below the plane of the figure. (Taken Dashed lines represent columns just above or below the plane 

from Ref. [^.) of the figure. (Taken from Ref g.) 
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3. Formation of vacancy /interstitial strings by slid- 



ing of polymers within columns in a columnar crystal of fi- 
nite-length polymers. 



FIG. 4. Formation of vacancy/interstitial strings by slid- 
ing of polymers within columns in a columnar crystal of fi- 
nite-length polymers. 
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Unlike dislocation lines, these strings (and loops) are not constrained to be planar: the lines can jump to any 
neighbouring lattice site as they traverse the crystal. Several horizontal jumps connecting a head to a tail are shown 
in Fig. ||. Note that left ward deflections of the interstitial segment connecting a head to a tail are accompanied by 
rightward deflections of the lines or columns themselves. A typical string can be approximated by an alternating 
sequence of straight segments and kinks joining the head of one column or polymer chain to the tail of another 
(see Fig. ||). Vacancy /interstitial strings are suppressed at low temperatures because they have a finite line tension, 
and hence an energy proportional to their length. At higher temperatures, heads and tails can move apart, forming 
variable- length strings that wander or "diffuse" perpendicular to their length by forming kinks. These strings thus 
resemble living polymers except that they are directed, on average, along the z-axis. In polymer crystals, the 
number of such strings is determined by the fixed concentration of heads and tails. In columnar discotic crystals, 
heads and tails can be created freely, and it is appropriate to treat their statistical mechanics in a grand canonical 
ensemble by introducing a head/tail fugacity, similar to the fugacity which controls defect concentrations in theories 
of vortex or dislocation unbinding transitions . We assume here that we can treat polymer crystals using the same 
formalism provided we tune the head/tail fugacity to achieve the fixed concentration determined by the mean polymer 
length. Long polymers imply a dilute distribution of heads and tails. We exclude, for simplicity, the possibility of 
hairpin excitations in polymer systems, which can be regarded as doubly quantized interstitial excitations leading to 
a higher energy. As we shall see, the sharp defect proliferation transition discussed in Ref. |^ is blurred when there 
is a finite concentration of heads and tails in equilibrium. 
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FIG. 5. Illustration of a vacancy string (thick dashed 
curve) joining a column head to another column's tail in a 
columnar crystal composed of long-chain polymers. 
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FIG. 6. Schematic of a defect string (composed of straight 
segments and kinks) wandering through the columnar crystal. 



Given an appropriate combination of parameters, namely, low line tension combined with head/tail and kink energies 
comparable to the temperature, the entropy of diffusion of the strings can overcome the line tension and lead to string 
proliferation, allowing heads and tails to separate to arbitrarily large distances. As in its bosonic counterpart, there 
exists off-diagonal long-range order in this phase, represented by 



lim (V'(r_L,z)i/'*(r_L',^')> 7^ 



(1.1) 



where tp and ip* are head and tail "destruction" and "creation" operators g, implying entanglement of lines on 
a macroscopic scale. If defects are absent or appear only in closed loops, the expression above would vanish as 
|rj^' — rj^l — > oo. Once defects proliferate, a line can wander to any other column and Eq. (1.1) has a finite limit. 
A crystal with proliferating vacancies and interstitials is an incommensurate phase — the magnitude of the smallest 
reciprocal vector G = in/^/^ao is no longer related to the areal density in the obvious way as p = \/3G^/8tt^ 
because the density differs from its defect- free value po = 2/V3ao ("^o being the lattice constant of the triangular 
lattice in cross-section). All crystals of pointlike atoms or molecules are trivially "incommensurate" in this sense 
— the corresponding pointlike vacancies and interstitials proliferate at any finite temperature. It is the anomalous 
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suppression of vacancies and interstitials and their organization into lines at low temperatures in columnar crystals 
which makes these materials unusual. 

In this paper we apply the physics of directed lines to vacancy/interstitial strings. With this in mind, we briefly 



review the elasticity theory of these systems in the next Section. In Sec. Ill we model a single string and estimate 



its transverse wandering. The form of this wandering is unchanged by coupling to phonon distortions of the lattice, 
as shown in Appendix |^ So is its magnitude, as calculated in Appendix In Sec. IV we apply the statistical 
mechanics of living polymers to an ensemble of directed strings and calculate their volume fraction, average length, 
etc. in the non-interacting limit. A simple quadratic-interaction model is presented in Section similar to the one 
discussed via the boson mapping in Ref . ||] , and we reproduce the results therein. Numerical calculations of the line 
tensions of various species of defects are presented in Sec. VI. The interaction potentials considered are repulsive and 
monotonic; we study simple power laws as well as a screened Debye-Hiickel interaction. We find many metastable 
species of vacancies. However, the lowest energy defect is always found to be the one with the highest symmetry in 
its category. For very short range interactions, this is the symmetric vacancy (Ve), whereas for most interactions the 
centered interstitial (/a) is most favoured. Appendix ^ contains details of the Ewald summation calculations for the 
potentials considered here. 



II. REVIEW OF ELASTICITY THEORY 



Before discussing defects in a columnar crystal, we review the aspects of elasticity theory common to all the systems 
mentioned in the Introduction. We consider lines or columns aligned along a common direction (z) up to thermal 
fluctuations, with crystalline order in any cross-section perpendicular to the columnar axis. In the case of flux lines, the 
average direction of alignment is imposed by an external field (H = Hz) and local deviations from this direction cost 
energy. With columnar crystals of long-chain molecules composed of covalently bonded nematogens or disk-shaped 
molecules cylindrically stacked via hydrogen bonds, or amphiphilic molecules in cylindrical micellar aggregates, the 
columnar axis represents spontaneously broken rotational symmetry. Therefore local deviations from the alignment 
direction are not penalized, but undulations of the column are. The rotational symmetry can, however, be broken 
by imposing an external field. In addition, the two-dimensional crystalline order resists shear and areal deformations 
perpendicular to the z-axis. 

Low-energy fluctuations of the system can be described by a "continuum" model that works for small ampli- 
tude, long- wavelength deformations [p|J^Jl0[|. The important fluctuations in this limit can be characterized by a 
two-dimensional displacement field u(rj^, z), representing the average deviation of lines in the {x, y) plane in a small 
region centered at (rj^, z). With it can be associated a local areal density change 5p/ = — Vj_ • u (po = 2/\/3ao) ^'^^ 
a local nematic director n = z + 1, with t = du/dz. The free energy of the system is a sum of nematic and crystalline 
contributions: 



^ — nematic "t" ^ crystal: 

To the lowest order in the fluctuations, these are given by 



J'nemat^c = ^ Jd^r [ifi(Vi • t)^ + K^iV ± X t)^ + K^id^t)^] 



and 



^ crystal — dz (Pr± 



pu,j^ + - A ( — 



(2.1) 



(2.2) 



(2.3) 



where Ki, K2, if 3 are the Frank constants for splay, twist and bend respectively, and A and p are the Lame coefficients. 
The matrix = {diUj -\- djUi)/2 is the linearized 2D strain field. In the presence of an external field iJz, one should 
add to T: 



J=ext = ^XaH^ Jdz Jd\^\t\ 



(2.4) 



where Xa is the anisotropic part of the susceptibility [|lO| . 

The last two contributions to J- are quadratic in the derivatives, and can be rewritten as 

^crystal + ^ext = 5 Jd^r [cii(V_L ' u)^ + C66(V_L X u)^ + d^id^u)^] 

+ fi (surface terms) 



(2.5) 
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where Cn = A + 2//, Cee = and C44 = XaH^P- The surface terms become important when there are defects within 
the bulk of the crystal, like vacancy/interstitial strings, represented by cuts joining column-end singularities in the 
field u(r_L , z) . Evaluating these terms over a cylindrical surface enclosing such a string yields the energy cost of the 
defect string: a line tension « /io^ due to the elastic distortion around the string, in addition to a core energy Ec 
per unit length (of the same order of magnitude) within the cylindrical core. 

Tnematic Can be further simplified if, as is often the case with ncmatic polymers, the splay and twist constants are 
small in comparison to the bend constant. Specifically, if Ki and K2 satisfy K i, 20^^ / \/ K^cn <C 1 0, then they can 
be neglected. For long- wavelength distortions along the columnar axis, the dominant free energy contribution is then 
K2,{dlu)^ in the absence of an external field. can be simply related to the persistence length Ip of the polymer as 
= ksTlpp. 

The statistical mechanics of defects in polymer liquid crystals has been discussd 
in detail by Selinger & Bruinsma lal. The presence of defects imposes a de- 
formation on the r = equilibrium configuration. In the case of a semi-infinite 
vacancy /interstitial string with a head or tail at the origin, this distortion fol- 
lows from minimization of the free energy above with respect to u(rj^, z) under 
the constraint 




Vj_ • u = ±Pq ^S{r±)9{z) + (non-singular terms) 



(2.6) 



FIG. 7. Distortion induced by 
a column end in the neighbouring 
columnar crystalline matrix. The 
distortion is confined to a vertical 
extent \z\ < a/XlTj" (shaded region) 
around the column end. 



where the ± sign refers to a column tail/head located at the origin. Since the 
planar distortion about a string has azimuthal symmetry in the continuum ap- 
proximation, Vj^ X u = 0. Hence, the only relevant terms in the free energy are 
the bend and bulk distortion terms (neglecting splay) . The resulting distortion 
around the column end spans a parabolic region about the radial direction (see 
Fig. 0) defined by 



(2.7) 



where Xl — ^/K^/cu is the length scale relating the distortions parallel and 
perpendicular to z. 

Selinger & Bruinsma also calculate the interaction energy between two column ends by superimposing the distortion 
created by each. They find the interesting result that a head and tail in a nematic medium attract weakly if they fall 
within each other's region of influence, as just described, but repel otherwise. However, in a columnar crystal (with 
non-zero shear modulus), the interaction is always a strong attractive linear potential due to the finite line tension 
associated with the string of distortions joining a head to a tail. 



III. WANDERING OF A SINGLE STRING 



Consider a single vacancy/interstitial string in a hexagonal columnar crystal of, say, polymer strands with lattice 
constant uq and monomer spacing c along the columnar axis z. For a discotic columnar liquid crystal, c is the spacing 
between oblate molecules along the column axis. For a flux line in a layered Type-II superconductor with magnetic field 
perpendicular to the layers, c is the layer spacing. If the string is vertical, the energy per unit length is of the order 
of /iOp (see Section ||) where /i is the in-plane shear modulus of the crystal. For a horizontal string, t± = £k/ clq where 
the kink energy Sk ^ K^^^fi^^^aQ k = K^/p being the bending rigidity. The ratio is tj_/tz ~ {k/ ii)^^^/ao ^ I* /a 
where /* is the kink size. Typically I* 3> oq, so that the strings are predominantly vertical, with few kinks. For flux 
lines on the other hand, the kink energy is g^^^p,^^^ao with g = C44/P, where C44 is the tilt modulus and p is the areal 
line density. The ratio is then {g/ p)^/"^ /ao- In highly anisotropic layered superconductors, this ratio can be small, 
favouring large, nearly horizontal defect excursions. We will for now work with nearly vertical strings, allowing for a 
gas of kinks sufficiently dilute so that the interaction between kinks can be ignored (see Fig. |). We thus assign to a 
string of vertical extent / and kinks an energy It + rikSk + 2eo where r = and eo is the energy of a polymer end. 
We expect that the results for defects with a high density of kinks would be qualitatively similar. 

In units such that fc^ = 1, the partition function of a string of length / is 

Zi-(l + ge-^-/^)'/'*e-'^/^ (3.1) 

where T is the temperature, and q is the two-dimensional co-ordination number of the lattice on which the defect 
string lives — for a symmetric vacancy this is the same as that of the original triangular lattice, q = 6, whereas for a 
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symmetric interstitial it is that of the dual honeycomb lattice, q = 3 (see Section VI). The above expression represents 
the freedom of the string to jump to any of the neighbouring lattice sites anywhere along its length. These transverse 
meanderings cause an entropic lowering of the free energy per unit length of the string: 



/i = lim -TlnZi// 

/ — ^oo 

^T-^e-'-l^ for e-^^-/^«l 
I* 



If iVfc is the total number of kinks, the average kink density is 

I ^ M + qe^^lT^ 

-'"Z^, for 



(3.2) 



I* 



< 1. 



(3.3) 



Thus, kinks are on the average Ik = l*e'^''^'^ /q monomers apart. The assumption of dilute kinks then translates into 
the condition l*nk <C 1, or, Sk » T, which can be rephrased as (|up)/a§ ^ 1 a condition clearly satisfied by a 

crystal below its Lindemann melting point. 

The above is a "diffusive" model for the string — if d denotes the horizontal end-to-end displacement, the mean 
square wandering is (|dp) = 2DI, where the "diffusion constant" D is given by 2D = UQUk- Consider a continuum 
description of the string in terms of a function rd{z), rd{z) being the transverse displacement. Provided the average 
slope \drd/dz\ is small, this "diffusive" wandering would correspond to an effective Hamiltonian of the form 



/ dz 
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dz 


+ T 


h 


2 





T 



(3.4) 



Here we have assumed that the string is wandering within a frozen crystal. However, the lattice around the 
vacancy/interstitial string responds to its presence by collapsing or expanding around it. For a straight string at 
r,^ — 0, the deformation u(r_L, z) is given by 



Ud(r_L, 



2tt r]_ 



(3.5) 



in the continuum description of the crystal, that is, away from the defect where the deformations are small. 51 is the 
area change due to the vacancy/interstitial, fl ~ Qq. The energy of this deformation has to be included in the energy 
cost of the defect string. Again invoking the continuum approximation, we assume that for a defect string with small 
average slope, the resulting deformation away from the string in any plane perpendicular to z would be approximately 
that resulting from a straight string at the location of the defect in that plane: 



u(r_L, z) ~ Ud(r_L - rd{z),z) 



(3.6) 



(In general u(rj^,z) would depend on the derivatives of Ydiz) as well.) Within this approximation, the distortion 
energy of the crystal with bending Frank's constant = Tlpp is, keeping terms up to fourth-order in the derivatives 
(see Appendix ^): 



Ai7i 
T 



h 



dz 



d'^Td 



dz^ 



dr„ 



dz 



(3.7) 



The se impar t an effective stiffness to the defect string and suppress transverse fluctuations over a length scale ^ 
aoy^ DK3/T ao\/lpnk- However, they do not change the long scale diffusive nature of the string. 

The lattice distortions renormalize the diffusion constant of the string when the symmetry direction of the crystal 
is externally imposed, as in the case of flux line s, or in a polymer crystal with an external field along the z-dircction. 
The tilt modulus C44 is then non-zero (Eq. (2.4)), and D is renormalized to Dp, where (see Appendix^ 
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For a dense vortex liquid this effect has been analyzed in detail by Marchetti ||13|] and D is found to be renormalized 
to a value independent of its bare value in the long- wavelength limit. The correction comes from convection of a 
tagged flux line along the local tangent-field direction. 

If a similar calculation is carried out for a crystal of spontaneously aligned long semi-flexible polymers (see Appendix 
^) , one finds a qualitatively different renormalization of D — the correction in the long- wavelength limit is proportional 
to its bare value, and 6D/D ^ 1.45(|up)/aQ ^ 3 % using cj^ ~ 1/50 [0 (cl is the Lindemann constant for melting of a 
columnar crystal). The correction is negligible. It can be ignored for another reason — the idea of convection of a line 
by the mean local field, although appropriate for a dense fluid, would not be applicable in a crystalline environment 
where diffusion can only occur through discrete jumps from column to column. Although thermal fluctuations are 
already implicit in the exponential factor in _D = coming from n^, defects in this case move only on a discrete 

lattice, without phonon fluctuations. 

To summarize this section, we characterize the statistical mechanics of a defect string with a head/tail energy So, a 
line tension r, and a diffusion constant D. The latter two can be combined in an effective chemical potential Jl = T/i^ 
per kink size (/*) of the string: 

= r(-r/T + uk) = qe-'"''^ - Efc/T, (3.9) 

with rik related to D through D — a^nk/l. Because nk is exponentially small, « —1*t/T « —l*^a^/T and is 
usually negative, which suppresses long vacancy & interstitial strings. Turning it positive would require raising the 
temperature and lowering the kink energy £fc, and is favored by a larger co-ordination number q. 

Although we have assumed a constant shear modulus, the presence of the defects themselves can drive it down 
exponentially with the defect concentration, as discussed by Carruzzo & Yu ||l5|. Thus, positive becomes possible 
when softening of the bare elastic constants with increasing defect concentration is taken into account. 



IV. STATISTICAL MECHANICS OF NON-INTERACTING STRINGS 



At any finite temperature, a crystal with a negative string line-chemical potential will contain a distribution of 
thermally excited vacancy and interstitial strings. Since the string energy is proportional to length in the non- 
interacting-kinks approximation, the equilibrium probability distribution would be an exponentially decaying function 
of length with mean determined by the line chemical potential, in the dilute string-gas limit where inter-string 
interactions can also be neglected [Q. In discotic crystals string heads and tails can be created as necessary. In a 
crystal of long polymers, the number of heads and tails is fixed by the mean polymer length. 

Let N be the total number of possible kink sites in the lattice, N = volume x p/l* , and Vi be 1/Nx the number 
of defect strings Winks long. Assuming that only one kind of defect string is present — those with the lowest line 
tension — we can write the defect free energy in terms of {Vi} as Q 

MiT^i}) = J2 ^^K2£o - IT fid) +Tj2NVi{lnri - 1) (4.1) 
I I 

Minimizing with respect to the {Vi} yields the expected exponential distribution: 

Vi = h^z^ (4.2) 

where z — e^"^, and the head/tail fugacity h = e~^°/'^ is expected to be small. For hexagonal columnar crystals of 
polymers, we work in a grand canonical ensemble and adjust eg so that the average head/tail concentration agrees with 
the fixed value determined by the mean polymer length. The head/tail concentration will be small if the polymers are 
long. For discotic crystals, the grand canonical ensemble is the natural one and the head/tail concentration fluctuates, 
with an average value determined by the fixed value oi h = e^^"^'^ , and the monomer fugacity z — e^'' < 1. The net 
defect volume fraction (t> is 



(1-^) 

The total number of strings Nd = Nus is given by the string density 



E'^' = '^'7r^- (4-3) 



E^'-T^ (4.4) 
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A defect monomer is most likely to be found in a string of mean length (in units of the kink size) 



lm = T^^ (4.5) 

The length distribution has an average at 21m, and a spread also of \/2lm- The form (|4.l| ) of the energy, linear in /, is 
really applicable only when Z 3> 1, so that end effects can be parametrized by the ^-independent constant Sq. Then, 
fid is close to 0, and the relation (/) ~ Uslmp holds. The asymptotic behaviours in the dilute and dense limits are as 
follows: 

4>={ (4.6) 



A string proliferation transition thus occurs at /z^ = in this model, corresponding to a temperature Td = rlk- In 
the limit Eq —> oo, it corresponds to the appearance of a supersolid phase |^] which is simultaneously crystalline and 
entangled, where infinitely long vacancy/interstitial strings facilitate the wandering and entanglement of lines in the 
crystalline phase. If the melting temperature Tm > Td, this supersolid/incommensurate solid phase will exist between 
Td and T^. 

The non-interacting approximation breaks down in the vicinity of Td as calculated here, and its estimate will have 
to be refined by including interactions. For finite eq, the sharp transition discussed in Ref. will be blurred, as 
discussed in Sec. 



-INTERACTION MODEL 



Interactions between polymer ends in a columnar crystal have been calculated by Selinger & Bruinsma pd| ] within 
the continuum approximation. Because of the uniaxial anisotropy, the interaction has a rather complicated form. 
The distortion due to an isolated head or tail plac ed a t the origin at in-plane distance r_L extends over a vertical 



extent \z\ ^ ^/XltZ where \l = \J K^/cn (see Eq. (2.7)). The resulting interaction between heads and tails falls as 
l/|zp for predominantly vertical separations z [\z\ 3> VAzTjT), and as for predominantly horizontal 

separations r_L. In polymer crystals, these contributions must be superimposed on the linear energy cost of the vacancy 
or interstitial string joining them. 

At low defect densities where the string length is much smaller than the average separation of string centers of 
mass, we have 1/4?-^^, i.e., \fXd\ ^ h?^^, and a string interacts with other strings as a head-tail dipole. The 

effective interaction between dipoles then falls off very rapidly, becoming short-ranged not only in the axial, but also 
in the radial direction. 

At the other extreme, the strings are long, which would happen in the vicinity of the head-tail unbinding transition 
and in the supersolid phase itself. End-interactions can then be neglected and the remaining interaction between 
effectively infinite strings becomes predominantly "radial" (i.e., perpendicular to z) provided the root mean square 
tilt with respect to the z axis is small. The defects are then non-interacting in the continuum model unless their 
anisotropy is taken into account. The interaction between defects with n-fold symmetry (n = 2, 3 or 6) falls off at 
least as fast as 1/r" (see Appendix ^). This interaction has an azimuthal dependence of the form cos nO or higher 
harmonics. The angular average vanishes, leading to an effective interaction which vanishes as an even higher power 
which is effectively short-ranged. As mentioned in the Introduction, the lowest-energy vacancy or interstitial defects 
for simple repulsive pair potentials in the radial direction are in fact of high (three-fold or six-fold) symmetry. 

We discuss here the simplest model for a short-ranged interaction — a repulsive cfp' model that has been treated 
earlier in Ref. |^ using a coherent state path integral representation which exploits an analogy with the quantum 
mechanics of two-dimensional bosons. The defect volume fraction (p corresponds to the mean square boson field 
amplitude (IV'P) in that description. Here, we reproduce the essential results witho ut re sorting to the sophisticated 



boson formalism. Upon adding a term m(/)^/2 to the free energy / = F/NT in Eq. (4.1) of the previous section, we 
find after minimization. 

As discussed in Ref. |^], the coupling u is an excluded volume parameter describing defect line repulsion. Thus and 
Nd have the same form as before, but with z replaced by an effective fugacity (: 



8 



z^C{z,(l)) = ze-'^'t', (5.2) 

so that 



^{hX)^h^j^^. (5.3) 



The volume fraction (/)(/i, z) now has to be solved for self-consistently from Eq. (5.3). Note that the effective chemical 
potential has been reduced by u(j) due to the repulsive interaction: 

/Xe// = In C = - (5.4) 
Accordingly, the mean string length l„i changes to 

l^ = -r^^—^. (5.5) 
In C U(p- fid 

The free energy of the distribution is / w —u(jP/2. 

The behaviour of the string volume fraction for h = and h ^ is illustrated schematically in Fig. ||. Four distinct 
regimes emerge, with the following asymptotic behaviours: 




This is again the dilute limit where heads and tails are tightly bound. 
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2. -1< < -(w/i2)i/3. 



1 



(5.7) 



These results are again identical to those for non-interacting strings. This correspondence is expected, because 
lA'dl > {uh'^Y/^ > u(j), therefore the effective chemical potential is still approximately /i^. The relation /id ~ 
~{uh?Y/^ marks the limit of validity of the non-interacting approximation, as we argued in the beginning of 
this section. As we approach this limit, we find for /i — + 0: (/), — > 0, whereas /,„ oo. Thus, the strings are 
still dilute, although lengthening. Note that the results in this regime coincide with those of Ref. ||] in the short 
and dilute strings limit. 



3. |/id| <C {uh'^Y^^ = He (^d around the transition which occurs for h = 0): 



1 + 1^ 
3 fie 



3 /ic. 



1 



3 Mc. 



(5.8) 



These results can be matched onto those in the non-interacting regime above by replacing figi with 



fJ-eff = -fJ-c + Aid/3 = -/ic 1 



(5.9) 



which is now dominated by the repulsive interaction: fieff ~ —ucj). The unphysical divergences of the non- 
interacting model have been suppressed and we find at the transition point: 



,,4/3 ■ 



,,1/3 ' 



1 



ul/3/^2/3- 



(5.10) 



Note that all quantities have interesting singularities in the limit h ^ 0. 

If the head/tail fugacity h is small, the defect volume fraction remains negligible at the transition, but the 
average string length grows large so that it could become greater than the inter-string separation, now given by 
1/^^/^. Indeed, l/cj)^^^ <^ Im H h \/u^ which would be true if polymer ends are highly unfavourable. 

This long & dilute regime interpolates between the short & dilute and the long & dense limits described in 
Ref. @. 

4. /id > ^JLc (Point B in Fig. |): 
In this limit, we have 



(5.11) 



The repulsion now keeps in check the string proliferation, and /ie// approaches as 1/ ^fjld. Thus, 



Md 

u 



1 



(5.12) 



This is the phase where strings are dense and entangled — is 0(1). These results also agree with Ref. [|. 

As the head/tail fugacity ft. 0, the intermediate regime 3 above (around /i = 0) shrinks to zero. At ft = 0, 
heads/tails are completely expelled, and we have a second-order phase transition at /i^ = with = for /i^ < 0, and 
growing as /i^ for /i^ > 0, as in Ref. This limit corresponds to the situation in thermally excited vortex lattices ||] 
because flux lines cannot start or stop within the sample. In the boson picture, ft acts like an external field coupled 
to the order parameter, injecting magnetic monopoles into the superconductor. 

We have neglected vacancy/interstitial loops, which exist even in the limit ft ^ 0. For finite ft, their contribution 
can be neglected near the transition because for long loops, the energy of a loop exceeds the energy of a string of 
the same vertical extent: Whereas a string of length / has energy iTintersUUai + 2eo (we expect interstitials to be 
the preferred defect at the transition in most cases), the energy of a vacancy-interstitial loop of the same length 
would be approximately / [ryacancy + nnterstiUai) ■ For large the difference Wyacancy - 2eo will strongly suppress 
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vacancy /interstitial loops. Because of this energetic barrier, loops cannot become arbitrarily large, and cannot cause 
entanglement over macroscopic scales. For /i = 0, as is the case for vortex matter, fluctuations in the low temperature 
phase are entirely in the form of loops and similar to vortex ring fluctuations in the Meissner phase. 

For systems with a finite axial length, the balance may be tilted in favour of long strings because the end penalty 
is removed if th e en ds move to the surface and the string threads the sample. For threading strings the expression for 



entropy in Eq. (O) is no longer valid because the freedom in the z-direction is lost. The remaining two-dimensional 



entropy can be ignored in a three-dimensional system, and we are left with 

/~-/irf0 + u0V2 (5.13) 
where </> now is also the areal fraction of defects; and one finds (j) ~ iJ-d/u, similar to region 4 discussed above. 



VI. NUMERICAL CALCULATION OF DEFECT LINE TENSIONS 



Line tension calculations require that we find the lowest energy lattice deformation associated with a vacancy or 
interstitial. These line tensions depend on the type of vacancy or interstitial, e.g., whether the defect sits in an 
environment which is two-, three- or six-fold symmetric. If thermal fluctuations out of this configuration are small 
enough to be described within a quadratic approximation, they decouple from the equilibrium configuration. Since 
these T — equilibrium defect configurations are composed of straight columns, the 3-dimensional deformation energy 
can be reduced to an effective 2-dimensional interaction energy V{r) per unit length between columns separated by 
distance r. The calculations can then be performed on a two-dimensional triangular lattice of points interacting with 
potential V{r). Thus, the defect energies in a two-dimensional Wigner crystal of electrons [ p^ would correspond 
to the line tensions of the corresponding string defects in a hexagonal columnar crystal of lines interacting with an 
effective radial l/r-potential per unit length. 

Such calculations have been carried out by several authors |^,|l6|,0. Whereas Refs. ||l^ and ||l^ have considered 
defects in a Wigner crystal of electrons {Vp{r) = 1/r), Frey et al. g have studied a modified Bessel-function potential 
Vfiir) = UQKQ{Kr) in the k ^ limit. Here k = A^^, where A is the Debye screening length in the case of long 
polyelectrolytes in an ionic solution, and the London penetration depth in the case of vortex lines in a type-II 
superconductor. The limit n corresponds to a long-range logarithmic interaction, whereas in the short-range 
limit Kao ^ 1 the interaction is exponentially decaying. Both Refs. Q and ||l^ dealt with long-range interactions (Inr 
and 1/r repectively), and found that the centered interstitial (see Fig. ^ has the lowest line tension. We denote the 
centered interstitial by CI, or by I3 when we want to stress its three-fold symmetry. The edge interstitial (denoted EI 
or I2) was found to be a saddle-point and buckled into a CI. The three- fold symmetric centered interstitial CI is the 
lowest energy interstitial defect over the entire range of interactions we studied. Among the vacancies, the two-fold 
symmetric crushed vacancy (denoted V2 or V2a — see Fig. ^) is the only stable one, the symmetric six-fold vacancy 
(Vg) being unstable to it. The long-range interactions between the energetically preferred types of interstitials and 
vacancies were found to be attractive for interstitials and repulsive for vacancies. 

To determine the correct type of microscopic defect to insert into the phenomenological considerations of Sees, 
p^-fv], we have extended the work of Frey et al. to the short-ranged regime of the i^Tp (Kr)-interaction, to which end 
we studied values of /tag from to 7 (7 being large enough to represent the short-range koq — *■ 00 limit) (Fig. |o|). 
The aim was to determine the point of cross-over from centered interstitials to vacancies as the lowest-energy defect, 
since it is known from simulations of short-range interactions (for a review, see Ref. fl^ ) that vacancies are preferred 
in this limit. In the same spirit, we have also extended the Coulomb interaction to power-law interactions l/r^ with 
exponent values ranging from p = (~ Inr) to p = 12 (Fig. |ll|). 

We checked our minimization procedure by first reproducing the results of Refs. and p7| for Inr and 1/r potentials 
respectively. As we move away from the long-range interaction limit na = 0, the metastable crushed vacancy (V2a) 
exchanges stability with the metastable split vacancy {SV), also of two- fold symmetry. Two metastable species, 
a three-fold symmetric vacancy (V3) and a two-fold symmetric vacancy (V2f,) crushed along the basis vector of a 
triangular unit cell, also exist, but are of higher energy. The differences in energy can be as small as one part in a few 
thousand. As the interaction gets shorter-ranged, V2& loses stability to V3 at kcq — 5.2, and the 3-fold deformation 
of V3 gets smaller so that it transforms continuously into Vq at kuq — 5.9. When Vq appears, the SV also loses 
stability to it. By the time I3 and Vq finally cross in energy, Vg is the only stable vacancy left. The crossing happens 
at surprisingly large parameter values, kuq ~ 6.9 for V^a (Fig. ^2|), and p ~ 5.9 for Vp (Fig. 13), each very close to 
the short-range limit. We thus find that the interstitial has a very wide range of stability, extending well into the 
short-ranged regime. 
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FIG. 9. Various defects obtained in a two-dimensional triangular lattice. The centered interstitial is the only stable interstitial 
defect. 
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FIG. 10. Defect energy as a function of the screening na 
for V{r) = Ko{Kr) at system size n = 4 {N ^ 480). Only the 
centered interstitial is shown, because the edge interstitial 
is always unstable to it. Various species of vacancies exist, 
within limited parameter ranges, very close in energy. Lines 
joining the data points are only an aid to the eye. 
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FIG. 11. Defect energy eis a function of the screening Ka for 
V{r) = at system size n = 5 (A'^ = 750). The apparent 
increase in energy with p (interaction getting shorter-ranged) 
would go away with proper normalization of the potential. 
Lines joining the data points are only an aid to the eye. 
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FIG. 12. Defect energies for V{r) — Ko{Kr), n = 4, on 
the log-scale, with respect to Vs/Ve, in order to illustrate the 
detailed structure of the energy diagram. The CI can be seen 
crossing Ve at na « 6.9. Lines joining the data points are only 
an aid to the eye. 
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FIG. 13. Defect energies for V{r) = l/r^, n = 5, on the 
log-scale, with respect to V^/Ve- The CI and Ve cross at 
p ~ 5.9. Lines joining the data points are only an aid to the 
eye. 



Following previous authors, the simulations were performed in an almost square (length-to- width ratio 5 : Sa/S) cell 
containing TV = 5n x 6n = 30n^ lattice points with n = 1 - 5 (rather than a more nearly square but bigger rectangle 
of, say, 7n x 8n (7 : 4-^/3) which would allow us to sample fewer number of system sizes n with a given computational 
limit on N). Fig. ^ corresponds to n = 3. 

A defect is introduced by adding or removing a particle, and then allowing the resulting configuration to relax. 
The difference between the energies of the relaxed defect configuration and the perfect lattice configuration gives the 
energy of the defect. There are two modifications to this simple calculation. We want the defect energy corresponding 
to the physical conditions of constant chemical potential or line density, so we rescale the cell dimensions (by changing 
the lattice constant uq) after inserting the defect to restore the system to its original density (following Ref. [|l7[). 
Moreover, since we would ideally like to study an infinite system, the large, but finite cell containing 30n^ particles 
is assumed to be repeated in all directions, so that we are effectively dealing with a periodic array of defects, or, an 
infinite lattice in the absence of a defect. The periodic boundary conditions maintain the average line density during 
the relaxation process. However, now the energy per cell also includes the energy of interaction of a defect with all 
its periodic images. As discussed earlier, this energy is finite, and by extrapolating its dependence on cell size n, i.e., 
inter-defect separation (« 5n), to large n, the energy of an isolated defect can be extracted P,p7[|. 

For short-ranged interactions, the energy calculation can be simplified. We introduce a cut-off interaction radius Vc 
where the interaction falls to a small fraction of its nearest-neighbour value. The interaction with the particles outside 
can be approximately accounted for by assuming a uniform density outside and integrating over it. The radius Tc is 
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chosen to make this correction smah compared to the total energy, say, less than 10~^ of it. Interactions within the 
shell are calculated explicitly. As long as Tc < i/2, L being the cell width, this short-range method should be very 
accurate. 

For long-ranged interactions such as Inr, l/r, or 1/r^, the above method breaks down, and we must resort to the 
Ewald summation technique [ p^j20| which yields an effective two-particle interaction that includes the interaction of 
one particle with all the periodic images of the other. This effective potential consists of a real space sum (corre- 
sponding to a screened interaction) and a reciprocal space sum (corresponding to the screening charge) . The division 
between the two is controlled by an Ewald parameter, and by a judicious choice of its value, the interaction can be 
made sufficiently short-ranged for both sums. We then employ cut-offs in both spaces, with values determined by the 
desired precision (see Appendix C for details). 
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.006835768 


6 


.005682901 


Ve 


Ve 


Ve 


Ve 


.006322377 


7 


.002788486 


V^e 


Ve 


Ve 


Ve 


.002771295 



TABLE I. Defect energies for V(r) = Ko{nr); ao — 1; system size n = 4 (TV = 480). The upper part corresponds to 
the Ewald Sum method for long-range interactions, the lower part to a simple cut-off method for short-range interactions. 
The centered interstitial and the symmetric vacancy cross at na ~ 6.9. Entries such as "V2a", "SV", "Va", "Ve" indicate an 
instability to a lower energy defect. 
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TABLE IL Defect energies for V(r) = 1/r^; ao = 1; system size n = 5 (A = 750). The Ewald sum technique was used to 



calculate the energies. The centered interstitial and the symmetric vacancy cross at p ~ 5.9. Entries such as "Va" and "Va" 
indicate an instability to a lower energy defect. 
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To find the minimum of the interaction energy as a function of the configuration of N particles, we use the conjugate- 
gradient method The forces are also needed for this method, and are easily derived from the energy and 
conveniently calculated along with it. 

The results for rt = 4 (480 particles) for V^a and for n — 5 for Vp (750 particles) are shown in Tables | and |l| and 
Figs. |l^ and |ll|. {n — 5 was computationally prohibitive for the long-ranged regime with kuq > 0). Note that, for the 
screened Bessel-function interaction, we find that calculations optimized for the long- and short-ranged regimes agree 
to within 1 part in 20, 000 at koq — 4. 

Moreover, we find that the interaction of a defect with all its periodic images is repulsive for defects with (even) 
two- and six- fold symmetry, and attractive for (odd) three-fold symmetry, consistent with Ref. As discussed in 
Refs. and [0, the true asymptotic form of the power law defect interaction probably isn't reached for the distance 
scales r ~ 20 — 30 lattice spacings studied here. 



VII. CONCLUSIONS 



We have studied factors contributing to the wandering of a vacancy or interstitial string defect in a hexagonal 
columnar crystal. A gas of such strings in the crystalline phase, interacting via short-range potentials, can proliferate 
via continuous or first-order transitions when the corresponding defect chemical potential changes sign, leading to 
a supersolid phase. The transition can be modified by the presence of vacancy-interstitial loops, especially in a 
system of finite thickness. We have also numerically calculated defect line tensions for two families of line interactions 
which interpolate between long- and short-ranged interaction potentials. In each case, we determine the point where 
interstitial and vacancy defects exchange stability. A complete accounting requires consideration of a variety of nearly 
degenerate vacancy configurations. At finite temperatures, the small energy differences between diffferent species will 
further lower the free energy of the vacancy through a gain in fluctuation entropy. The interstitial itself can fluctuate 
between the centered and edge configurations. The point where vacancies and interstitials exchange stability will shift 
at finite temperatures due to entropic effects of this kind. In the context of long-range potential calculations, we show 
in an Appendix how to extend the Ewald summation to the modified Bessel function potential Kq{x). 
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APPENDIX A: CALCULATION OF ENERGY OF DISTORTION DUE TO A DEFECT STRING 



As described in Ref. minimization of the free energy (2.1) with the constraint (2.6) yields the following equation 
for u(rj_, z): 

A|5> - (Vi • u) = -y±S{rj_ - Vd) (Al) 
Po 

Td being the in-plane location of the defect string (assumed straight for now). Upon assuming a solution of the form 
u = — ^Vj_'0, we have the scalar equation 

i-Xldt + Vi)V = <5(rj. - Vd) (A2) 

For the straight string, the solution is 

■(/'(r_L,2) oc ln|r_L - r<j|, or, u(rj_, z) cx -— (A3) 

Irj- - rd\ 

with proportionality constant ~ Oq. 

Now consider a wandering string with a dilute concentration of ki nks, described on average by rd{z) (see Fig. |^). 
Upon inserting this z-dependence into the right hand side of Eq. ( [A2| ), we see that the resulting i/j inherits the 
fiuctuations of rd{z). If Iz represents the smallest wavelength in rd{z), the two terms on the LHS of Eq. ( |A2[ ) compare 
as Xj^/l^ vs. I/cq, or as I* vs. Iz where I* — VAlOq is of the order of the kink length. Since the meandering of the 
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defect string occurs on a length scale much larger than the kink size, the first term should be negligible compared to 
the second, and we can set 



i{}{r±,z) oc ln|r_L -rd(z)| 



(A4) 



as a reasonable approximation. 

The elastic energy of a defect of length L can now be written as 



E, 



defect 



TzL + Ek 



dz 
ao 



dYd{z) 



dz 



d^rK3 



(A5) 



representing contributions from line tension, kinks, and the bending energy of the distorted crystal (zero for a straight 
string). The primed integral here excludes the core of the string: a region of radius ao around it. It can easily 
be evaluated for u(r, z) = Ud{r± — rii{z),z) and reduces to the form in Eq. (3.7), accurate up to fourth order in 
the derivatives. The second term, on the other hand, leads to the term {g/2) Jdz\drd/dz\'^ in Eq. (3.4). For long 
wavelengths, the additional contribution from the third term is irrelevant in comparison, being of higher order in the 
derivatives. The length scale at which it becomes important is obtained by balancing the two terms: K^/l^ ^ T/D 
or, U ~ ^/K^D/T. 



APPENDIX B: RENORMALIZATION OF D BY DEFECT-PHONON COUPLING 




FIG. 14. Illustration of the coupling 
between a defect string and the lattice dis- 
tortion. In this case, the change in the po- 
sition of the vacancy string (thick dashed 
curve) is equal and opposite to the change 
in the phonon displacement field. 



In Section |l| we described the wandering of a defect line along the z- 
axis by a "diffusio n" c onstant D = aQnk/2, corresponding to an effective 
Hamiltonian (Eq. (3^)) Hdefect — (5/2) Jdz\drd/dz\'^ , rd{z) describing the 
in-plane position of the defect string, with g = T/D. To incorporate the 
effect of lattice fluctuations on the diffusion of the defect string, we modify 
Hdefect to 



Hdefect-phonon — 2 



dVd 

±^+t(r„.) 



(Bl) 



where the expression in brackets now represents the deviation of the va- 
cancy/interstitial string with respect to the local director 



du 

dz 



(B2) 



Fig. ^ illustrates the case of a vacancy string, which we shall assume for 
the remainder of this appendix. 



It is easy to derive the diffusion equation for the partition function Z(rc(, Tq; z, 0) corresponding to the above 
Hamiltonian (the {u(rj^, z)}-dependence in Z has been omitted for convenience): 



d^Z - (t • Vj_)Z 



L»V^ Z 



(B3) 



Z represents the probability density for the defect position; — t is the "convective velocity" for this density. It can 
also be thought of as an (imaginary) vector potential acting on a particle of mass g in two dimensions, with z the 
time-like coordinate. 

Defining the propagator G{r±,z) — Z{r±, z)9{z), d(z) being the step function, G obeys 

(9, - DWl)G{r^,z) = d'^^\r^)S{z) + 1 • V^G (B4) 
The bare propagator Go corresponds to ignoring the convective influence of the medium. Thus, Go satisfies 

{d,-D\7l)Goir^,z)^S<^'\r) (B5) 
Fourier-transforming rj^ — > k (space- like) and z lu (time- like). 
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The renormalized diffusion coefficient Dji will be calculated from the average of G over the phonon degrees of freedom 
using the definition 



in the limit |k|,a; 0. Upon denoting k = (k, w), Eq. (B4) becomes 

Go\k)G{k) = 1 + [ ik'-t(fc-fc')G(fc') 

Jk' 

The symbol denotes / cPk/{2'K)^. Eq. ( p^ ) can be expanded in a perturbation series: 
G{k) = Goik) + Go{k) [ ik' • t(fc - fc')Go(fc') 



+ Go{k) / ik' ■t{k-k')Goik') / ik"-t(fc'-/c")Go(fc") + 

Jfc' Jk" 

To calculate the thermal averages of products of t = —iuju, we need 

K(fc)) - 0, 

Mk)u(,{k')) = [SL{k)Vtp{k) + ST{k)Vlp{k)\ 5^^\k - k') 
= S^p{k)5^^\k - k') 
where the correlation functions parallel (L) and perpendicular (T) to k are 

T 



SL/rik) 



and the projection operators are 'P^fj{k) — kakp/k'^ and V^pik) — Sa/s — T-'apik). Therefore 



(G(fc)) =Go(fc)-Go(fc) 
Diagrammatically, this series is represented in Fig. |15|. 
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FIG. 15. Diagrammatic representation of the series expansion of the propagator for the defect probability density. The 
average is over the phonon degrees of freedom. 
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All diagrams of type (b), (c) reducible to the one-loo p dia gram (b) can be summed to give (G(fc)) "'^ « Go(fc) -'^ + 
V{k) where V{k) is the term in square brackets in Eq. (B13). Then, Dj^ — D + SDji where 

5Dn^\im.UjlV{k) (B14) 

With the assumption Xi^/tD^o^)'^ ^ Ij in other words, kink size I* ^ Ik, the kink separation, 6Dji can be evaluated 
to give 

where we have imposed a cutoff by replacing the hexagonal Brillouin zone by a circle of radius A = which has 



the same area. If cge cn, the fluctuations are mostly transverse, and 

5Dr ^ (|uH 47r 
D ^ al 5V3 



(B16) 



Since ^^"i ^ ^ c\ where the Lindemann ratio for hexagonal columnar crystalline lattices is empirically known to be 
4 ~ 1/50 [||, we find ^ < 3%. 



APPENDIX C: INTERACTION BETWEEN DEFECTS WITH A-FOLD SYMMETRY 

We assume the defect string to be straight, so that only the planar elastic deformation energy, T crystal (Eq. is 
relevant. In the continuum model discussed, this energy is isotropic in the strains Uij. The stress Uij — SJ-'crystai/^Uij 
can be expressed in terms of a biharmonic stress function x [ p2| : V^V^x = 0, as ct^ — tik^jidkdix [f-rj is the 
two-dimensional anti-symmetric tensor, ei2 — 1). 

For a dislocation with Burger's vector b, x = ~Kh x rj^ Inrj^ where K ~ fi{X + ii)/tt{X + 2fi) in terms of the Lame 
coefficients. 

We construct x for ri-fold symmetric vacancy/interstitial by treating it as a superposition of n dislocations 
(bound) symmetrically placed a distance d ^ ao apart (such that the volume of the defect is f2 ~ ndao), with Burger's 
vectors separated by 27r/n in orientation. The resulting stress function has a form satisfying 

Vi,.™(2..»fc,--l|„.^) (CI, 

where are coefhcients of 0(1). 

The interaction energy of two such defects, located at ri and r2 respectively, can be written in terms of their stress 
functions as 

Z^i2(ri2 = r2 - ri) = ^ y d\i_Vlx{ri_ - ri)Vix(r± - v^). (C2) 

For ri2/d ^ 1, the leading term in the interaction comes from the convolution of the (5- function with the fc = 1 term 
(in other words, this is the cost of the volume change produced by one defect in the stress field of the other), therefore 
it is of the form cos n9 12 /r 12^'' . 

Specifically, we find for vacancies (the sign is reversed for interstitials) : 



n 


relative orientation 


U{rj_) in units oi Kfl/d'' 


2 
2 
3 
3 
6 


parallel 
perpendicular 

parallel 
anti-parallel 


-cos 2«+cos4« 

-2cos2e+cos4e 

— 3 cos 66 

— 2cos3y 1 3 cos 6a 
(r/d)^ ' (r/d)^ 
-2cos6d 
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APPENDIX D: THE EWALD SUM FOR V{R) = KoiKR) AND 1/7?^ 



Let (/)(r) be the two-body interaction potential between charges qi in the given system so that the total interaction 
energy is 

The simulated system consists of N particles in a cell repeated to generate an infinite system. Then, the energy per 
cell can be written as 

1 ^ 

U = Uo^ + ^Y.Yl ^^j^a ) (D2) 

where 

n G 



lini [va{r) 

r— >0 



0(r)] 



(D3b) 



Here the sum over n consists of all real space lattice vectors of the lattice generated by a cell of area A, and the sum 
over G goes over the corresponding reciprocal lattice vectors. The lattice is rectangular (almost square) in our case, 
which makes it easy to list these vectors, ipa is the long-range part of the interaction 0, so that (/> — ?/;q, is a screened, 
short-ranged interaction, ipa is the Fourier transform. The amount of screening is controlled by the Ewald parameter 
a. If we take a to be large enough so that the real space sum can be truncated at Tc = L/2 within the desired 
precision, then we can drop the sum over n p3t . However, this means including more short-range components 
in the screening charge distribution, so that it spreads to higher reciprocal vectors. The cutoff in reciprocal space is 
again determined by the precision required. 

Since we shall be considering N particles, each with qi — 1, we have J^i = ^ ^'i^d liY — ■ Rearranging 
the sums in U and noting that N / A = p which is constant throughout the calculation, we can rewrite U as 



U = Uref + U„ 



where mm 



w 1 

U^nt « E E ['^('■^j ) - V-aCr.,)] + ^ E 



^cosG-r.j +1^^ 



sin G • r , 



(D4) 



(D5a) 



U, 



ref 



N N ~ 

-— lim V^„(r) + — p lim V„(k) 



(D5b) 



Note that Uref is explicitly proportional to N in this form. This is the form we use in our calculations. For interactions 
whose long-range integral diverges (such as l/r^ with p < 2), a uniformly spread background of equal and opposite 
charge is assumed, so that the second term in Uref should contain 



lim 



^„(k)-0(k) . 

We can now proceed to the special potentials we are interested in. For the power-law potential 

0(r) ^l/r'P^ 

we take il 



1 

r(p72) Jo 



dt 



(D6) 



(D7) 
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^ r(p/2)i„. rJ' r(p/2) ^ ^ 



so that the screened interaction is 



1 7 (p/2, (gr)^) 

(^^-^-^'^^^^^F r(p/2) ■ 



(r and 7 are complementary incomplete Gamma functions p5[|.) 
The Fourier transform is (d = 2): 



*,i^)-^(it'r(i^Jl) 1 (mo, 



r(p/2) Vfc/ V 2 ' \2a 

Also, 



lim ipaik) 



p-d r(p/2) 



p> d 



(D12) 



For p < d, the above expression corresponds to limfc^o i^aik) — (j){k) . For p = d, both forms would diverge, however, 

they would be independent of a, and since the defect energy is a difference of energies, this term would cancel out. 
The force on particle j, fj = —\7jU, can also be written as a sum of real space and reciprocal space contributions: 



(D13a) 



where 



- 2 

^3 



r(p/2 + l,(ar,,f 



(D13b) 



is the sum of forces on particle j due to particle i and all its images, and 

= 1 E C^) (Y. co*^ ^ • G • r, - sin G • r, j 



cos G • r 



G 



(D13c) 



is the sum of forces on particle j due to all images of itself. 

The co-ordinates r here are normalized such that oo = 1 . When we change A'' to Nd = iV ± 1 and rescale uq to 
a — \J Nd/N after inserting a defect, we chose to keep r normalized with respect to a, so that it picks up a factor of 
a. If we also scale a by 1/a, the product ar remains unchanged (as does G • r), so that we can keep using the original 
values of a and G in Uint & f, and scale the result by in the end. In Uref we have to use the scaled value of a 
along with Nd, and subtract the energy of the perfect lattice scaled by Nd/N. 

On the other hand, if we do not scale a, Uref cancels out, but a has to be replaced by aa in Uint & f ■ 

For the modified Bessel function interaction 



0(r) = Koinr), (D14) 

(where k represents na because of the normalization of ?■), the k — case, corresponding to (p{r) ^ — Inr, has been 
treated by Frey et al. |Q. To extend this to k > 0, we perform an expansion similar to that of Mokross & Silva ]2^ ] 
for a Yukawa potential. Writing the potential in integral form, 

0(r) = K^inr) = - / -e''e'^ (D15) 
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We choose the screened interaction to be 

-i'a{r) = - / — e « 



1 



2 7i s 

Expanding the exponential within the integral in a Taylor series about k = 0, we get 

" (-1)" / K N2n 



(D16) 



(D17) 



0, only the first term, Ei ((ar)^), is non-zero. For k > we have an alternating series, and its convergence 



where En{x) is the exponential integral function {En+i{x) — x'^T{~n, x)). 
For K 

has to be taken into account in determining the optimum value of a (in addition to the required precision and the cell 
size). For large values of k, not only are a large number of terms needed in this series to reach the desired precision, 
the optimum value of a is also large due to convergence rquirements, increasing the cutoff in reciprocal space, so that 
the computation time increases dramatically. We were able to carry these calculations to k = 4, where it matched 
the results from the short-range method to 1 part in 20, 000. 
We also need the following quantities (in d = 2): 



(D18) 



27r /- K 
lim ipa{k) ^e^l.^ 



(D19) 



At K = we take 



Also, 



lim lim 



V'a(fc)-0(fc) 



TT 

'2^ 



lim Va(r-) « -7/2 + lna-ln(K/2) +0(k^) 



(D20) 



(D21) 



(D22) 



The In k, term cancels in the defect energy, so that the limit k — > is again well-defined. Similarly, if we did not 
subtract limfc^o 4'{k), we would have an extra term 27r/K ^, which too would cancel out. 

The expression for the force is similar to Eqs. ( D13| ) where f,;^ is a series similar to 0(r) — ij^air), with each 

((ar)^) replaced by 2a^£'„ ((ar)^) r. On scaling, ^a{k) has to be recalculated because it does not simply scale 
as a power law. 

The special functions used here were all calculated to an accuracy of 10^^^ according to routines taken from 
Ref. pl| ] . Since the power-law calculations were mostly carried out on integral values of p, the gamma functions were 
only needed for integral or half- integral orders, in which case certain recursion relations can be used [^. We used 
the fastest method for each order. 
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